/**
 * @file ebex_histogram.h
 *
 * @date Sep 26, 2011
 * @author seth
 *
 * @brief This file is part of EBEX State Estimator, created for the EBEX project
 *
 * This software is copyright (C) 2011 Columbia University
 *
 * EBEX State Estimator is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * EBEX State Estimator is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with EBEX State Estimator; if not, write to the Free Software Foundation, Inc.,
 * 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

#ifndef EBEX_HISTOGRAM_H_
#define EBEX_HISTOGRAM_H_

#include <stdint.h>
#include <stdbool.h>



typedef struct emxArray_real_T
{
	double *data;
	int32_t *size;
	int32_t allocatedSize;
	int32_t numDimensions;
	bool canFreeData;
} emxArray_real_T;

/*
 *
 */
void histc(const emxArray_real_T *X, const emxArray_real_T *edges, emxArray_real_T *N, emxArray_real_T *BIN)
{
	bool b0;
	int32_t nbins;
	int32_t xind;
	uint32_t outsize[2];
	int32_t bin;
	uint32_t uv0[2];
	bool eval_single = FALSE;
	uint32_t m;
	int32_t exitg1;
	uint32_t low_i;
	uint32_t low_ip1;
	uint32_t high_i;
	uint32_t mid_i;

	nbins = edges->size[0];

	outsize[0] = 1U;
	outsize[1] = (uint32_t) nbins;

	xind = N->size[0] * N->size[1];
	N->size[0] = 1;
	xind = N->size[0] * N->size[1];

	N->size[1] = (int32_t) outsize[1];
	bin = (int32_t) outsize[1] - 1;

	for (xind = 0; xind <= bin; xind++)
	{
		N->data[xind] = 0.0;
	}
	for (xind = 0; xind < 2; xind++)
	{
		uv0[xind] = (uint32_t) X->size[xind];
	}
	xind = BIN->size[0] * BIN->size[1];
	BIN->size[0] = 1;
	xind = BIN->size[0] * BIN->size[1];
	BIN->size[1] = (int32_t) uv0[1];

	e_memset(BIN->data, 0, sizeof(double) * uv0[1]);

	eval_single = FALSE;
	if (nbins > 1)
	{
		m = 2;
		do
		{
			if (m <= (uint32_t) nbins)
			{
				if (edges->data[m - 1] < edges->data[m - 2])
				{
					xind = N->size[0] * N->size[1];
					N->size[0] = 1;
					xind = N->size[0] * N->size[1];
					N->size[1] = (int32_t) outsize[1];
					bin = (int32_t) outsize[1] - 1;
					for (xind = 0; xind <= bin; xind++)
					{
						N->data[xind] = rtNaN;
					}
					xind = BIN->size[0] * BIN->size[1];
					BIN->size[0] = 1;
					xind = BIN->size[0] * BIN->size[1];
					BIN->size[1] = (int32_t) uv0[1];
					bin = (int32_t) uv0[1] - 1;
					for (xind = 0; xind <= bin; xind++)
					{
						BIN->data[xind] = rtNaN;
					}
					break;
				}
				else
				{
					m++;
				}
			}
			else
			{
				eval_single = TRUE;
				break;
			}
		} while (1);
	}
	else
	{
		eval_single = TRUE;
	}
	if (eval_single)
	{
		nbins = X->size[1];
		xind = 0;
		for (m = 1U; m <= (uint32_t) nbins; m++)
		{
			bin = 0;
			if ((!(edges->size[0] == 0)) && (!isnan(X->data[xind])))
			{
				if ((X->data[xind] >= edges->data[0]) && (X->data[xind] < edges->data[edges->size[0] - 1]))
				{
					low_i = 1U;
					low_ip1 = 2U;
					high_i = (uint32_t) edges->size[0];
					while (high_i > low_ip1)
					{
						mid_i = (low_i >> 1U) + (high_i >> 1U);
						if (((low_i & 1U) == 1U) && ((high_i & 1U) == 1U))
						{
							mid_i++;
						}
						if (X->data[xind]
								>= edges->data[mid_i - 1])
						{
							low_i = mid_i;
							low_ip1 = mid_i + 1U;
						}
						else
						{
							high_i = mid_i;
						}
					}
					bin = (int32_t) low_i;
				}
				if (X->data[xind] == edges->data[edges->size[0] - 1])
				{
					bin = edges->size[0];
				}
			}
			if (bin > 0)
			{
				N->data[bin - 1] += 1.0;
			}
			BIN->data[xind] = (double) bin;
			xind++;
		}
	}
}

#endif /* EBEX_HISTOGRAM_H_ */
